#|(defmeth homals-proto :compute (
               &key
               (data (send self :data))           ; data matrix
               (ndim (send self :p))              ; # computing dimensions
               (active-variables                  ; list of active variables
                 (send self :active-homals-variables))
               (active-categories                 ; list of active categories
                 (send self :active-categories))  ; for the active variables
               (object-labels                     ; labels of all objects
                   (send self :object-labels))
               (variable-labels                   ; labels of all variables
                   (send self :variable-labels))
               (category-labels                   ; labels of all categories
                   (send self :category-labels))
               (eps-0 (send self :eps-0))         ; function change
               (eps-1 (send self :eps-1))         ; solution change
               (itmax 200)                        ; # iterations
               (speed t)                          ; Speed or Memory
               x                                  ; Z for bootstrapping
               bootstrap                          ; Performing Bootstrapping?
               bootnum                            ; Bootstrap Iteration Number
               (original-data 't))                ; Using original data?

(let* ((n (send self :n))
       (e (select data (iseq n) active-variables))
       (m (length active-variables))
       (fit-0 0)
       (fit-1 0)
       (itel 0)
       (row-ind-list (if (and bootstrap (not original-data))
                         (send self :sample-row-ind-list)
                         (send self :row-ind-list)))
       (col-ind-list (if (and bootstrap (not original-data))
                         (send self :sample-col-ind-list)
                         (send self :col-ind-list)))
       (dlist (if (and bootstrap (not original-data))
                  (send self :dlist-samp)
                  (send self :dlist)))
       (w (if (and bootstrap (not original-data))
              (send self :weights-samp)
              (send self :weights)))
       (x (if x x (make-random-orthonormal w n ndim m))))
  (loop
    (let* ((y (make-cat-quant x col-ind-list dlist ndim))
           (new-fit-0 (* n ndim))
           (new-fit-1 0)
           (z (make-obj-scores y row-ind-list w ndim)))

     (setf new-fit-0 (- new-fit-0 (/ (sum (diagonal (matmult (transpose y) 
                        (apply #'bind-rows (* dlist (row-list y)))))) m)))
     (setf z (* (sqrt n) (gram-schmidt (homals-center z w) (/ w m))))
     (setf new-fit-1 (/ (sum (combine (matmult (transpose z)
                                      (apply #'bind-rows (* w (row-list x))))))
                        (* m n ndim)))

     (if (or (and (> eps-0 (abs (- new-fit-0 fit-0)))
                  (> eps-1 (abs (- new-fit-1 fit-1))))
             (> itel itmax))
         (progn
               (send self :d-m nil)
               (let ((c (make-array (list ndim ndim) :initial-element 0))
                     (y (make-cat-quant z col-ind-list dlist ndim))
                     (k-j-sumlist (cons 0 (cumsum (send self :k-j-list))))
                    )
                 (mapcar #'(lambda (a b) 
                          (let* ((ymat (select y (iseq a (1- b)) (iseq ndim)))
                                 (dmat (diagonal
                                           (select dlist (iseq a (1- b)))))
                                 (d (/ (matmult (transpose ymat) dmat ymat) n)))
                            (setf c (+ c d))))
                     (butlast k-j-sumlist) (rest k-j-sumlist))
                 (let* ((f (elt (eigen c) 1))
                        (k (apply #'bind-columns f))
                        (z (matmult z k))
                        (y (matmult y k)))
                  (send self :z z)
                  (send self :y y)
                 (mapcar #'(lambda (a b)
                          (let* ((ymat (select y (iseq a (1- b)) (iseq ndim)))
                                 (dmat (diagonal
                                           (select dlist (iseq a (1- b)))))
                                 (d (/ (matmult (transpose ymat) dmat ymat) n)))
                            (send self :d-m d)))
                     (butlast k-j-sumlist) (rest k-j-sumlist)))

                 (if (and bootstrap (not original-data))
                   (format t
                       "\nBootstrap Iteration ~3d Complete...\n\n" bootnum)
                   (format t "\nDone....\n\n"))
                 (return)))
       (progn
         (format t "~%Iteration ~4d, Loss Measure ~,10f, Change measure ~,10f"
                 (1+ itel) new-fit-0 new-fit-1)
         (setf fit-0 new-fit-0)
         (setf fit-1 new-fit-1)
         (setf itel (1+ itel))
         (setf x z)))))))|# ;removed PV


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;  Set selection lists
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmeth homals-proto :set-inds ()
 (let* ((data (send self :data-matrix))
        (n (send self :n))
        (active-variables (send self :active-homals-variables))
        (active-categories (send self :active-categories))
        (e (select data (iseq n) active-variables))
        (m (length active-variables))
         row-ind-list
         col-ind-list)
      (send self :k-j-list nil)
      (dotimes (j m)
         (let* ((g (make-indicator (elt (column-list e) j)
                                       (elt active-categories j)))
               (row-inds (mapcar #'(lambda (x) (which (= 1 x))) (row-list g)))
               (col-inds (mapcar #'(lambda (x) (which (= 1 x)))
                            (column-list g)))
               (kj (send self :k-j-list)))
           (setf kj (if kj (sum kj) 0))

           (if (= j 0) 
               (progn 
                 (setf row-ind-list row-inds)
                 (setf col-ind-list col-inds))
               (progn 
                 (setf row-ind-list (mapcar #'(lambda (x y)
                       (if y (append x (+ kj y)) x))
                            row-ind-list row-inds))
                 (setf col-ind-list (append col-ind-list col-inds))))
           (send self :k-j-list (array-dimension g 1))))
      (send self :row-ind-list row-ind-list)
      (send self :col-ind-list col-ind-list)
      (send self :dlist (mapcar #'length col-ind-list))
      (send self :weights (mapcar #'length row-ind-list))))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Output routines
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmeth homals-proto :homals-numerical-output (output-file)
  (let (
        (data (send self :data-matrix))
        (z (send self :z))
        (y (send self :y))
        (d-m (send self :d-m))
        (av (send self :active-homals-variables))
        (ac (send self :active-categories))
        (ol (send self :object-labels))
        (vl (send self :variable-labels))
        (cl (send self :category-labels))
        (p (send self :p))
        (n (send self :n))
        (w (send self :weights))
        (m (send self :m)))
    (with-open-file
     (output output-file :direction :output)
     (format output "~a~%" (make-string (+ 27 (* 13 p)) :initial-element #\=))
     (format output "Object Scores~%")
     (format output "~a~%" (make-string (+ 27 (* 13 p)) :initial-element #\=))
     (dotimes (i n)
         (format output "~20a *** A " (elt ol i))
       (dotimes (s p)
         (format output "~12,8f " (aref z i s)))
       (format output "~%"))
     (format output "~a~%~a~%" (make-string (+ 27 (* 13 p)) 
                        :initial-element #\=) #\C-L)
     (format output "~a~%" (make-string (+ 36 (* 13 p)) :initial-element #\=))
     (format output "Category Quantifications~%")
     (format output "~a~%~%" (make-string (+ 36 (* 13 p)) :initial-element #\=))
     (format output "~a~%" (make-string (+ 36 (* 13 p)) :initial-element #\-))
     (dotimes (j m)
      (format output "~a~%" (make-string (+ 36 (* 13 p)) :initial-element #\-))
      (if (find j av)
           (progn
           (format output "Variable ~20a      (ACTIVE)~%" (elt vl j))
           (format output "~a~%" (make-string (+ 36 (* 13 p)) 
                           :initial-element #\-))
           (let* ((r (elt (column-list data) j))
                  (u (sort-data (remove-duplicates r :test 'equal)))
                  (g (make-indicator r u))
                  (y (make-category-quantifications z g))
                  (d (column-sums g))
                  (k (length d)))
             (dotimes (l k)
               (if (find (elt u l) (elt ac (position j av)) :test 'equal)
                   (format output "~20a *** ~4d *** A "
                           (elt (elt cl j) l) (elt d l))
                   (format output "~20a *** ~4d *** P "
                           (elt (elt cl j) l) (elt d l)))
               (dotimes (s p)
                 (format output "~12,8f " (aref y l s)))
               (format output "~%"))))

           (progn
            (format output "Variable ~20a     (PASSIVE)~%" (elt vl j))))
       (format output "~a~%" (make-string (+ 36 (* 13 p)) 
                      :initial-element #\-)))

     (format output "~a~%" #\C-L)
     (format output "~a~%" (make-string (max 43 (* 13 p)) 
                    :initial-element #\=))
     (format output "Discrimination Measures~%")
     (format output "~a~%~%" (make-string (max 43 (* 13 p)) 
                    :initial-element #\=))

     (let ((e (make-array (list p p) :initial-element 0)))
       (dotimes (j m)
         (format output "~a~%" (make-string (max 43 (* 13 p)) 
                    :initial-element #\-))

       (if (find j av)
           (progn
             (format output "Variable ~20a      (ACTIVE)~%~%" (elt vl j))
             (let ((d (elt d-m (position j av))))
               (setf e (+ e d))
               (dotimes (s p)
                        (dotimes (u p)
                                 (format output "~12,8f " (aref d s u)))
                        (format output "~%"))))
           (progn
               (format output "Variable ~20a     (PASSIVE)~%" (elt vl j))))
         (format output "~a~%" (make-string (max 43 (* 13 p)) 
                         :initial-element #\-)))
       (setf e (/ e (length av)))
       (format output "~a~%" (make-string (max 43 (* 13 p)) 
                   :initial-element #\-))
       (format output "Average Discrimination Measure~%")   
       (format output "~a~%" (make-string (max 43 (* 13 p)) 
                   :initial-element #\-))
       (dotimes (s p)
                (dotimes (u p)
                       (format output "~12,8f " (aref e s u)))
              (format output "~%"))
       (format output "~a~%" (make-string (max 43 (* 13 p)) 
                   :initial-element #\-))))))



(defmeth homals-proto :make-opt-dat-mat (output-file)
  (let (opt-mat 
        (data (send self :data-matrix))
        (z (send self :z))
        (n (send self :n))
        (m (send self :m)))
   (dotimes (i m)
      (let* ((r (elt (column-list data) i))
             (u (sort-data (remove-duplicates r :test 'equal)))
             (g (make-indicator r u))
             (y (bind-columns (first 
                   (column-list (make-category-quantifications z g)))))
             (mati (matmult g y)))
        (if opt-mat (setf opt-mat (bind-columns opt-mat mati))
                    (setf opt-mat mati))))
   (with-open-file (f output-file :direction :output)
     (let ((opt-rows (row-list opt-mat)))
       (dotimes (i n)
         (mapcar #'(lambda (x) (format f "~6,4f " x))
              (coerce (elt opt-rows i) 'list))
         (format f "~%"))))))


(defmeth homals-proto :write-object-scores (output-file)
  (let ((z (send self :z)))
    (with-open-file (f output-file :direction :output)
       (dotimes (i (array-dimension z 0))
          (dotimes (j (array-dimension z 1))
              (format f "~7,5f " (aref z i j)))
          (format f "~%")))))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Computing methods for an ALS iteration.
;;
;; It would be useful to have versions of this in C,
;; and load them dynamically.
;;
;; The same thing is true for the Gram-Schmidt and the Procrustus
;; modules below, because that is where most of the computing takes
;; place.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun make-category-quantifications (x g)
  (let* ((d (geninv (column-sums g))))
    (matmult (diagonal d) (matmult (transpose g) x))))

(defun make-discrimination-measures (y g)
  (let* ((d (column-sums g)))
    (/ (matmult (transpose y) (matmult (diagonal d) y))
       (array-dimension g 0))))



(defun make-cat-quant (x col-ind-list dj p)
 (apply #'bind-rows (mapcar #'(lambda (a b)
               (/ (column-sums (select x a (iseq p))) b)) col-ind-list dj)))

(defun make-obj-scores (y row-ind-list w p)
 (apply #'bind-rows (mapcar #'(lambda (a b)
               (/ (column-sums (select y a (iseq p))) b)) row-ind-list w)))

(defmeth homals-proto :set-sample ()
 (let* ((n (send self :n))
        (samp (sample (iseq n) n t)))
    (send self :get-boot-rows samp)
    (send self :get-boot-columns samp)
    (send self :get-boot-dlist)
    (send self :get-boot-weights)))


(defmeth homals-proto :get-boot-rows (samp)
  (send self :sample-row-ind-list (select (send self :row-ind-list) samp)))

(defmeth homals-proto :get-boot-columns (samp)
 (let* ((n (send self :n))
        (numcat (sum (send self :k-j-list)))
        (sample-col-ind-list (repeat nil numcat))
        (sample-row-ind-list (send self :sample-row-ind-list)))
   (mapcar #'(lambda (x y) (mapcar #'(lambda (z)
                (setf (elt sample-col-ind-list z)
                      (cons y (elt sample-col-ind-list z)))) x))
      sample-row-ind-list (iseq n))
   (send self :sample-col-ind-list sample-col-ind-list)))

(defmeth homals-proto :get-boot-dlist ()
  (send self :dlist-samp (mapcar #'length (send self :sample-col-ind-list))))

(defmeth homals-proto :get-boot-weights ()
  (send self :weights-samp (mapcar #'length (send self :sample-row-ind-list))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Utilities
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun column-sums (x)
  (mapcar #'sum (column-list x)))

(defun column-means (x)
  (mapcar #'mean (column-list x)))

(defun row-sums (x)
  (mapcar #'sum (row-list x)))


(defun make-indicator (x &optional (z nil set))
  "Args: sequence
Elements of SEQUENCE are either numbers or strings.
Returns a dummy with sorted category values."     
(let* (
      (y (if set z
           (sort-data (remove-duplicates x :test 'equal))))
      (m (length y))
      (n (length x))
      (a (make-array (list n m) :initial-element 1))
      (b (make-array (list n m) :initial-element 0))
      )
  (if-else (outer-product x y #'equal) a b)
))

(defun marginals (x)
  "Args: sequence
SEQUENCE is a sequence of numbers or strings. Different entries are
sorted and counted."
  (mapcar #'sum (column-list (make-indicator x)))
  )

(defun number-of-values (x)
  "Args: sequence
Elements of SEQUENCE are either numbers or strings.
Returns the number of different values."
  (length (remove-duplicates x :test 'equal))
  )

(defun make-random-orthonormal (w n ndim m)
  "Args: w n ndim m
Makes a matrix of order N x M with standard normals,
then centers and w-orthogonalizes it."
  (let ((z (make-array (list n ndim) :displaced-to
                       (coerce (normal-rand (* n ndim)) 'vector))))
    (* (sqrt n) (gram-schmidt (homals-center z w) (/ w m)))))

(defun q-r-decomp (x w)
"Args: X W 
X is decomposed as QR, with Q w-orthonormal."
 (let ((v (geninv (sqrt w)))
       (y (first (qr-decomp x 't))))
   (apply #'bind-rows (* v (row-list y)))))


(defun gram-schmidt (x w)
  "Args: X
X is decomposed as KS, with K w-orthonormal and S upper-triangular,
returns K."
  (let ((y (chol-decomp (matmult (transpose x) 
              (apply #'bind-rows (* w (row-list x)))))))
    (matmult x (inverse (transpose (first y))))
    ))

(defun homals-center (x w)
  "Args: X W
X is a matrix and W is a list of weights.  Returned is a matrix Z such that
u'WZ=0"
  (let (
        (n (first (array-dimensions x)))
        (mu (/ (matmult w x) (sum w)))
        )
    (- x (outer-product (repeat 0 n) mu #'+))
    ))

(defun geninv (x &optional (ozo 1))
  (if (compound-data-p x)
      (map-elements #'geninv x)
    (if (= 0 x) ozo (/ x))))


(provide "core.lsp")
